# Load CIFTI data file
cii <- read_cifti(mscfile, drop_data = FALSE, trans_data = T)
# Make brainstructure index
### Notes: The brainstrucure index lets us filter out anatomical structures based on an
### index (this mirrors the cifti packages in MATLAB)
cii$brainstructureindex <- as.matrix(NA, dim(cii$data)[1])
for(i in 1:length(cii$BrainModel)){
startindx <- as.numeric(attributes(cii$BrainModel[[i]])$IndexOffset + 1)
endindx <- as.numeric(attributes(cii$BrainModel[[i]])$IndexOffset +
attributes(cii$BrainModel[[i]])$IndexCount)
cii$brainstructureindex[startindx:endindx] <- i
}
# Load parcel and community assignment of each vertex (provided with MSC data)
### Notes: Community is used for identifying parcels into community
### (i.e., sub-network like Default Mode Network)
parcel <- read_cifti(parcel_file)
parcel <- as.matrix(parcel$data)
comm <- read_cifti(comm_file)
comm <- as.matrix(comm$data)
# ==== Check dimension of cifti data (volume/frame x vertices)
dim(cii$data) # ~ 64k vertices, includes subcortical volumes
## [1] 65890 818 1
dim(parcel) # surface only, excluded medial wall
## [1] 59412 1
dim(comm) # surface only, excluded medial wall
## [1] 59412 1
# What are the labeled brain structures in the cii file?
### Notes: CIFTI data contains the surface (cortx left and right) and subcortical
### structures based on volumetric data. The labels should contain the left &
### right coritcal surface, and individual subcortical labels.
cifti_brain_structs(cii)
## [1] "CIFTI_STRUCTURE_CORTEX_LEFT"
## [2] "CIFTI_STRUCTURE_CORTEX_RIGHT"
## [3] "CIFTI_STRUCTURE_ACCUMBENS_LEFT"
## [4] "CIFTI_STRUCTURE_ACCUMBENS_RIGHT"
## [5] "CIFTI_STRUCTURE_AMYGDALA_LEFT"
## [6] "CIFTI_STRUCTURE_AMYGDALA_RIGHT"
## [7] "CIFTI_STRUCTURE_CAUDATE_LEFT"
## [8] "CIFTI_STRUCTURE_CAUDATE_RIGHT"
## [9] "CIFTI_STRUCTURE_CEREBELLUM_LEFT"
## [10] "CIFTI_STRUCTURE_CEREBELLUM_RIGHT"
## [11] "CIFTI_STRUCTURE_HIPPOCAMPUS_LEFT"
## [12] "CIFTI_STRUCTURE_HIPPOCAMPUS_RIGHT"
## [13] "CIFTI_STRUCTURE_PALLIDUM_LEFT"
## [14] "CIFTI_STRUCTURE_PALLIDUM_RIGHT"
## [15] "CIFTI_STRUCTURE_PUTAMEN_LEFT"
## [16] "CIFTI_STRUCTURE_PUTAMEN_RIGHT"
## [17] "CIFTI_STRUCTURE_THALAMUS_LEFT"
## [18] "CIFTI_STRUCTURE_THALAMUS_RIGHT"
# ==== Take the Left & Right cortex only
### Notes: subcortical data are not sorted into the community assignments provided by MSC
### data, so only the cortical surface is being analyzed.
cdata <- as.matrix(cii$data[cii$brainstructureindex==1 | cii$brainstructureindex==2,,])
# cdata <- cii$data[,,]
dim(cdata)
## [1] 59412 818
u_parcel <- unique(parcel)
u_parcel <- u_parcel[u_parcel!=0] # Remove parcel 0 and order parcel by their number
# ==== Mask out bad frames/volumes from data (i.e., scrubbing; Power et al. 2012, 2014)
### Notes: The time-mask (tmask) that identifies bad frames are provided with MSC data.
tmask <- read.table(tmask_file)$V1
ctmask <- cdata[,as.logical(tmask)]
# ==== Extract mean time series from each parcel itno a matrix (Parcel x Frame)
tp <- matrix(0, length(u_parcel), sum(tmask))
tp_notmask <- matrix(0, length(u_parcel), length(tmask))
for(i in 1:length(u_parcel)){
tp[i,]<- colMeans(ctmask[which(parcel==u_parcel[i]),])
tp_notmask[i,]<- colMeans(cdata[which(parcel==u_parcel[i]),])
}
# ==== Order matrix by parcel number
tp <- tp[order(u_parcel),]
tp_notmask <- tp_notmask[order(u_parcel),]
superheat::superheat(tp,
heat.lim = c(-20, 20),
heat.pal = c("black","white"),
grid.hline = FALSE,
grid.vline = FALSE,
title="Mean Time series of each parcel (bad frame removed)")
superheat::superheat(tp_notmask,
heat.lim = c(-20, 20),
heat.pal = c("black","white"),
grid.hline = FALSE,
grid.vline = FALSE,
title="Mean Time series of each parcel (bad frames retained)")
r <- cor(t(tp)) # Correlation matrix between all nodes
z <- psych::fisherz(r) # Fisher's z-transform: 0.5 * log((1+r)/(1-r))
diag(z) <- 0 # Set diagonal to '0'; not informative
superheat::superheat(z,
y.axis.reverse = TRUE, # Used to make origin (0,0) on top left corner
heat.lim = c(-.2, .4),
heat.pal = rev(brewer.rdylbu(100)),
heat.pal.values = c(0, 0.15, 0.25, 0.75,1),
grid.hline = FALSE,
grid.vline = FALSE,
title="Node x Node Correlation Matrix (z-transformed)")
parlabel <- data.frame(parcel_num=sort(u_parcel), community=NA, comm_label=NA, comm_shortlabel=NA)
plotlabel <- read.csv("../data/systemlabel_MSC.txt", header=F,
col.names = c("community","comm_label","color","comm_shortlabel"))
for(i in 1:length(u_parcel)){
comm_i <- unique(comm[which(parcel==sort(u_parcel)[i])])
parlabel$community[i] <- comm_i
parlabel$comm_label[i] <- as.character(plotlabel$comm_label[is.element(plotlabel$community, comm_i)])
parlabel$comm_shortlabel[i] <- as.character(plotlabel$comm_shortlabel[is.element(plotlabel$community, comm_i)])
}
superheat::superheat(X = z,
y.axis.reverse = TRUE,
membership.rows = parlabel$comm_shortlabel,
membership.cols = parlabel$community,
left.label.col=plotlabel$color,
bottom.label.col=plotlabel$color,
extreme.values.na = FALSE,
heat.lim = c(-.2, .4),
heat.pal = rev(brewer.rdylbu(100)),
heat.pal.values = c(0, 0.15, 0.25, 0.75,1),
title="Parcel x Parcel Correlation Matrix (z-transformed)")
# ==== Setup positive matrix plot
z_pos <- z
z_pos[z<0] <- 0
ss_pos <- superheat::superheat(X = z_pos,
y.axis.reverse = TRUE,
membership.rows = parlabel$comm_shortlabel,
membership.cols = parlabel$community,
left.label.col=plotlabel$color,
bottom.label.col=plotlabel$color,
extreme.values.na = FALSE,
heat.lim = c(0, .3),
heat.pal = parula(20),
heat.pal.values = c(0, 0.5, 1),
title="Node x Node Positive Correlation Matrix (z-transformed")
# ==== Setup negative matrix plot
z_neg <- z
z_neg[z>0] <- 0
ss_neg <- superheat::superheat(X = z_neg,
y.axis.reverse = TRUE,
membership.rows = parlabel$comm_shortlabel,
membership.cols = parlabel$community,
left.label.col=plotlabel$color,
bottom.label.col=plotlabel$color,
extreme.values.na = FALSE,
heat.lim = c(-.3, 0),
heat.pal = rev(parula(20)),
heat.pal.values = c(0, 0.5, 1),
title="Node x Node Negative Correlation Matrix (z-transformed")
gridExtra::grid.arrange(ggplotify::as.grob(ss_pos$plot), ggplotify::as.grob(ss_neg$plot),
nrow=1)
ss_smooth_pos <- superheat::superheat(X = z_pos, smooth.heat = T, smooth.heat.type = "mean",
y.axis.reverse = TRUE,
membership.rows = parlabel$comm_shortlabel,
membership.cols = parlabel$community,
left.label.col=plotlabel$color,
bottom.label.col=plotlabel$color,
extreme.values.na = FALSE,
heat.lim = c(0, .3),
heat.pal = parula(20),
heat.pal.values = c(0, 0.5, 1),
title="Node x Node Positive Correlation Matrix (z-transformed")
ss_smooth_neg <- superheat::superheat(X = z_neg, smooth.heat = T, smooth.heat.type = "mean",
y.axis.reverse = TRUE,
membership.rows = parlabel$comm_shortlabel,
membership.cols = parlabel$community,
left.label.col=plotlabel$color,
bottom.label.col=plotlabel$color,
extreme.values.na = FALSE,
heat.lim = c(-.3, 0),
heat.pal = rev(parula(20)),
heat.pal.values = c(0, 0.5, 1),
title="Node x Node Negative Correlation Matrix (z-transformed")
gridExtra::grid.arrange(ggplotify::as.grob(ss_smooth_pos$plot), ggplotify::as.grob(ss_smooth_neg$plot),
nrow=1)
## Threshold matrix to 4%
z4 <- z_pos
z4[z < quantile(z, 0.96)] <- 0
net <- graph.adjacency(adjmatrix = z4, mode = "undirected", diag = F, weighted = T)
parlabel$id <- 1:nrow(parlabel)
parlabel$color <- NA
u_comm <- unique(parlabel$community)
for(i in u_comm){
parlabel$color[parlabel$community==i] <- as.character(plotlabel$color[plotlabel$community==i])
}
V(net)$id <- parlabel$id
V(net)$community <- parlabel$community
net <- simplify(net, remove.multiple = F, remove.loops = T)
pnet <- plot(net, layout=layout_with_fr, vertex.label=NA, vertex.size=5,
vertex.color=parlabel$color, alpha=.6)
if (!require("NetworkToolbox", character.only=T, quietly=T)) {
devtools::install_github("AlexChristensen/NetworkToolbox")
}
##
## Attaching package: 'NetworkToolbox'
## The following objects are masked from 'package:igraph':
##
## betweenness, closeness, degree, diversity, strength,
## transitivity
## The following object is masked from 'package:dplyr':
##
## desc
library(NetworkToolbox)
p <- participation(A = z4, comm = parlabel$community)
# Each node's PC calculated using positive & negative edge
# Negative edges were taken out in previous steps, so PC caculated with all-edges and positive-edges are the same.
head(p$overall)
## [1] 0.8375559 0.3965046 0.6669420 0.5441093 0.6031861 0.3480727
head(p$positive)
## [1] 0.8375559 0.3965046 0.6669420 0.5441093 0.6031861 0.3480727
# PC based on negative edges should all be zero (not usable in this case).
head(p$negative)
## [1] 0 0 0 0 0 0
# Coloring nodes based on PC.
gcol <- grey.colors(n=nrow(z))
plot(net, layout=layout_with_fr, vertex.label=NA, vertex.size=5,
vertex.color=gcol[order(p$positive)], alpha=.6)
parlabel$pc_4td <- p$positive
ggplot(parlabel, aes(x=pc_4td)) +
geom_histogram(bins=20) +
xlab("Participation Coefficient (4%)") +
ggtitle("Participation Coefficient (4%) across entire network") +
theme_bw()
ggplot(parlabel, aes(x=pc_4td)) +
facet_wrap(~comm_shortlabel) +
geom_histogram(bins = 20) +
xlab("Participation Coefficient (4%)") +
ggtitle("Participation Coefficient (4%) across each sub-network") +
theme_bw()
Use a density plot to visualize distributions where there are small number of nodes.
ggplot(parlabel, aes(x=pc_4td)) +
facet_wrap(~comm_shortlabel) +
geom_density() +
xlab("Participation Coefficient (4%)") +
ggtitle("Participation Coefficient (4%) across each sub-network") +
theme_bw()
parlabel$degree_4td <- degree(A = z4) # calculate degree for each node at 4% edge density
ggplot(parlabel, aes(x=degree_4td)) +
facet_wrap(~comm_shortlabel) +
geom_histogram(bins = 20) +
xlab("Degree (4%)") +
ggtitle("Degree (4%) across each sub-network") +
theme_bw()